Fluctuation dissipation ratio in an aging Lennard-Jones glass 
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By using extensive Molecular Dynamics simulations, we 
have determined the violation of the fluctuation-dissipation 
theorem in a Lennard-Jones liquid quenched to low temper- 
atures. For this we have calculated X{C), the ratio between 
a one particle time-correlation function C and the associated 
response function. Our results are best fitted by assuming 
that X{C) is a discontinuous, piecewise constant function. 
This is similar to what is found in spin systems with one step 
replica symmetry breaking. This strengthen the conjecture of 
a similarity between the phase space structure of structural 
glasses and such spin systems. 

PACS numbers: 61.43. Fs, 61.20.Lc , 02.70.Ns, 64.70. Pf 



Obtaining information on the phase space structure of 
glassy systems is a very difficult challenge. By definition, 
relaxation times in a glass are so long as to preclude 
equilibration within an experimental (or numerical) time 
scale, except perhaps for very small systems (l]j^. Ex- 
ploration of phase space in these systems is necessarily 
incomplete, and the results from any experimental inves- 
tigation cannot be expected to be representative of a well 
defined statistical ensemble. Hence, although many con- 
jectures have been formulated concerning the structure 
of phase space in glassy systems Q , very little is actually 
known. 

A promising route, that might to some extent bypass 
this intrinsic difficulty, is the idea that relevant informa- 
tion on phase space structure is encoded in the nonequi- 
librium dynamics of glassy systems. This idea was ac- 
tively developed in the field of spin glasses but its 
extension to the field of structural glasses is more recent 
[0-|lo|. Among the important quantities that can be in- 
vestigated in a nonequilibrium system is the so called 
fiuctuation-dissipation ratio X, defined in the following 
way. Consider an observable A whose normalized auto- 
correlation function will be denoted by C, and let R be 
the response function of A to a field H conjugate to A. 
Then, for a system in equilibrium at temperature T, C(t) 
is related to the response function R{t) of the system to 
H by the usual fiuctuation dissipation theorem (FDT), 
R{t) — In a system that is out of equilibrium 
(e.g., a system that has been quenched at t = to a 
low temperature) the property of time translation invari- 
ance is lost, and C and R become functions of two time 



variables, e.g. C{t',t) = {A{t')A{t)). A formal way of 
generalizing the usual FDT consists in writing, for t' > t 
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Rit',t) = -—x{t\t) 



d C{t',t) 
'di 



(1) 



which in this form is merely a definition of X. The impor- 
tance of this "FDT violation factor" X{t',t) was recog- 
nized in the context of mean field theories of spin glasses 
]ll[ , where it appears that X{t' , t) has the properties dis- 
cussed below. For this discussion it is useful to consider 
the situation in which the system is driven out of equilib- 
rium at time t = Q, then aged for a waiting time tw after 
which the measurement of the time correlation functions 
C{tu, -I- r, t^,) are started. For mean field models exhibit- 
ing glassy behavior it has been shown that in the limit 
of long times, i^,, t — > oo, X{tyj + r, tyj) is a function of 
the correlation function C only, i.e. 



X{t^ -l-r, tyj) = x{C{t^ +T,t^)), 



(2) 



where x is now a function of one variable. In this asymp- 
totic limit two regimes can be distinguished. If is 
kept fixed, C{tw + T,tw) eventually becomes a function 
of r only. The limiting value of this function for r — > oo 
is the Edwards- Anderson parameter associated with ob- 
servable A, which is usually denoted by qea- Obviously 
qea vanishes for an equilibrium system, and differs from 
zero in a nonergodic system. For 1 > C > Qea, we have 
x(C) = 1, meaning that the FDT holds. This means 
that for time differences that are small compared to the 
waiting time t^, the response of the system is similar 
to that of an equilibrium system, in spite of the fact 
that only a restricted part of phase space is explored. 
Nonequilibrium, or "aging" features, show up in a differ- 
ent limit, namely for t > t^. In this limit, the correlation 
function depends on both t^, and r in a nontrivial way, 
typically C(iu, -|-t, t^,) = F{h{tw + T)/h{T)), where h{x) 
is a monotonically increasing function. In this "aging" 
regime, C < Qea, and x{C) < 1. The system starts to 
sample a larger portion of phase space, but this sampling 
is an out of equilibrium process, and does not obey the 
equilibrium FDT. 

An important property of x{C), again discovered in 
the context of mean-field spin glasses, is that the gen- 
eral structure of this function is identical to that of the 
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function Xstat (9) obtained by inverting the integral of the 
Parisi function qstat (x) ■ The latter reflects the probabil- 
ity distribution of overlaps between replicas of the same 
system, and does not involve any dynamical considera- 
tion At present, the similarity between these two 
functions has escaped physical understanding, although 
a formal justification has recently been proposed [ p^ . 
This similarity between the two functions is, neverthe- 
less, believed to be a general feature. If this is the case, 
the implication is that the study of a dynamical quan- 
tity such as x{C) provides indirect information on the 
structure of phase space. So far x(C) was determined for 
the 3d Edwards- Anderson model |14|, for ferromagnetic 
coarsening fl^ j, for p-spin models in 3 dimensions [^ , 
for the 3d Ising spin glass and a string in a disordered 
medium fl^ ] , confirming in each case the general features 
of mean-field predictions. In this Letter, we show that 
an accurate determination of x{C) for a structural glass 
model is indeed possible, using standard simulation tech- 
niques. 

The system we study is a 80:20 mixture of 1000 
Lennard-Jones particles, with interaction parameters 
that prevent crystallization In the following, we 

shall use as length, energy and time units the standard 
Lennard-Jones units aAA (particle diameter), caa (inter- 
action energy), and t = (mAC^^/48eAA)^^^, where ttia 
is the particle mass and the subscript A refers to the 
majority species [Q. The system has been described in 
detail elsewhere, and its equilibrium (high temperature) 
properties have been fully characterized. At the reduced 
density p = 1.2, a "computer glass transition" is found 
in the vicinity of T = 0.435 and the slowing down of the 
dynamics seems to be described well by Mode-Coupling 
theory |jl^. A first study of the aging behavior of the 
correlation functions at low temperatures has also been 
published recently 0]. 

In order to obtain a fluctuation dissipation ratio, we 
need to compute C and R for the same observable. Pre- 
vious work 01 focused on the aging behavior of the inco- 
herent scattering function for wavevector k: 

Ck{t^+T,t^) ^ l^e*-(--^(*»+^)-^^(*"». (3) 
j 

In order to compute the associated response function, 
we use the following numerical approach. A fictive 
"charge" e = ±1 is assigned randomly to each parti- 
cle. An additional term of the form J^j Cj^(i"j)> where 
V{r) = Vq cos(k-r) is a small {Vq < ksT) external poten- 
tial, is then added to the Hamiltonian. It is then easy to 
show that, if one averages over several realizations of the 
random charge distribution, the time-correlation function 
of the observable Ak — J2j exp(ik • rj{t)) is the inco- 
herent scattering function. The procedure to generate 
the response function associated to Ck is thus straightfor- 
ward: For a given realisation of the random charge distri- 



bution, the system is equilibrated at a high temperature 
(T = 5.0), and quenched at t = to the desired final 
temperature Tf. The evolution is followed with the field 
off until a waiting time t^, then the field is switched on 
and the response Ak{tw -l-r, iiu) is monitored. The same 
procedure is repeated for several (7 to 10) realisations of 
the charge distribution, in order to get the response func- 
tion. The quantity we obtain by this procedure is then 
an integrated response function M{tw + t, i™), defined 
as: 

{Ak {tw + T, tm 

)) = Vo / R{t^+T,t)dt (4) 



VoM{t^, + T, t^). 



(5) 



This procedure was carried out for three different val- 
ues of the final temperature Tf, namely Tf = OA,Tf — 
0.3 and Tf = 0.1. The amplitude of the external poten- 
tial was chosen in such a way that a linear response is ob- 
tained at each temperature. For Tf = 0.4, Vq — 0.2 while 
for Tf = 0.1 Vb = 0.05. The wavevector was k = 7.25, 
the location of the main peak in the structure factor. The 
runs had a length of 5 • 10^ time steps, corresponding to 
100,000 time units. 

Typical data for the integrated response and for the 
correlation function is shown in Fig. p, for Tf = 0.4 and 
k = 7.25. The way by which the correlation function 
and the integrated response are related to each other can 
be understood well by means of a parametric plot of M 
versus C, as shown in Fig. ^ for different values of t^ and 
two different Tf . If the generalized fluctuation theorem 
holds, it is easily checked that M can be written as a 
function of C, with 



M(C) = -■ 
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x{c)dc 



(6) 



ksT jc 

From Fig. ||, it is clearly seen that the usual FDT with 
X = 1 is very well verified for short times, i.e. values of C 
close to 1 , in that the curves are linear and have slope — 1 . 
For longer times a break in the curves is observed, i.e. the 
FDT is violated. Some transient effects are perceptible 
for the shorter waiting time, but they tend to disappear 
with increasing tw This violation is compatible with 
the Ansatz (^, since the parametric curves obtained for 
different waiting times superimpose satisfactorily. 

For the regime in which the FDT is violated, the M{C) 
curve can have different forms |^,^: Domain growth 
models predict that M{C) is a constant, whereas mean 
field models predict a linear dependence, for the case of 
"one step" replica symmetry breaking, and a more gen- 
eral dependence for the case of continous replica symme- 
try breaking. 

As can be seen from Fig. ^ a good fit to the result- 
ing M versus C curve is obtained with a piecewise lin- 
ear function, which corresponds to a piecewise constant 
xiC): 
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c(C) = 1 for C > qb, x{C) = m < 1 for C < 96, 



(7) 



where is the value of C at which the mentioned break 
in the curves is observed. Such a dependence has, e.g., 
been found for mean-field "p-spin" models [|. Thus our 
results give support to the hypothesis first formulated by 
Kirkpatrick and coworkers and revived by Parisi , 
that structural glasses belong to the same "universality 
class" as mean-field p-spin models. 

From Fig. |^a we can read off m ~ 0.62 and ~ 0.6. 
The latter value is clearly smaller than the plateau value 
for the correlation function in Fig. |l|, qEA — 0.8. This 
means that the FDT appears to hold even for times at 
which the system is no longer time translationally invari- 
ant, a feature which is not predicted by current theories 
of aging. 

Finally, the dependence of x on the final temperature 
can be investigated. To explore this dependence we have 
also done simulations at Tf = 0.3 and Tf = 0.1. In all 
cases, we find that the M vs. C plot can be fitted well by 
two straight lines. Our results for m as a function of T 
are thus given by: Tf = 0.4, m = 0.62 ± 0.05; Tf = 0.3, 
m = 0.45 ± 0.05; Tf = 0.1, m = 0.2 ± 0.1 Within the ac- 
curacy of our data these values of m are compatible with 
a linear dependence on Tf, quite similar to that found 
by Parisi for a soft-sphere system. Such a linear de- 
pendance (m(T/) ^ Tf) corresponds to a constant "fluc- 
tuation dissipation effective temperature" Te// — Tf/m. 
The later concept, introduced in could help ratio- 
nalize the older "Active temperature" idea. 

We mention that, in his analysis of the fluctuation 
dissipation relation, taking as an observable the mean 
squared displacement, Parisi found that rn{Tf) can be 
approximated by m(Tf) = Tf/Tc for Tf < T^, where T^ 
is the "mode coupling critical temperature" of the sys- 
tem under study. In our case, Tc ~ 0.435 [|l^, so that our 
data is in contradiction with such a simple dependence of 
m on T . Our results are much more similar to the ones 
found by Alvarez et al. for the p-spin model in that these 
authors found for a temperature a bit below Tc a value 
of m which is significantly smaller than 1. The reason 
for this difference might be related to the much smaller 
waiting times used in Ref. |9| [0. In any case, it is not 
clear why the mode coupling critical temperature should 
play a particular role in the present analysis. If the same 
type of simulations would be carried out at a tempera- 
ture slightly above Tc, we expect that interrupted aging 
will be observed, so that at short tyj a violation of FDT 
occurs. As increases, equilibrium will progressively be 
approached, and the M versus C plot will approach a 
straight line with slope —1. Hence the main difference 
between T > and T < T^ will be that for T > T^ 
the M versus C plot depends on tw, as has been shown 
for the Sherrington-Kirkpatrick model in three dimen- 
sions iQ. However, for T close to, but above, Tc this 
dependence will be so weak that it can be neglected for 



all practical purposes. In terms of the "effective temper- 
ature" Teff = Tf/m, our system falls out of equilibrium 
above Tc, so that we can expect Tg// to be larger than 
Tc - which is indeed the case. 

In summary, we have shown that the fluctuation dis- 
sipation ratio of a supercooled liquid out of equilibrium 
can be computed with good accuracy from MD simula- 
tions. Several nontrivial features predicted by the theory 
of mean-field spin glasses, beginning with the existence 
of a waiting time independent function x{C), seem to be 
present also in the model structural glass we study. Our 
data is compatible with a stepwise constant x{C), which 
would correspond to a phase space structure similar to 
that of spin systems undergoing one step replica sym- 
metry breaking. This means that phase space is divided 
by high barriers into different valleys each of which has 
the same statistical properties. (The case of continuous 
replica symmetry breaking corresponds to a case that the 
valleys are organized in a hierachical way.). In any case, 
finding a nonzero value of m seems to be a clear indi- 
cation that a "domain growth" picture is not applicable 
to our model. A quantitative comparison between theo- 
retical predictions and simulation results, similar to the 
work carried out for testing mode coupling theory, would 
be required in order to fully clarify the situation with 
respect to nonequilibrium dynamics. However, unfortu- 
nately such calculations are currently not yet feasible. 
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FIG. 1. Correlation function C{tn, + T,tn:) (dashed lines) 

and integrated response function M(tw + T,t^) (solid lines) 
for Tf = 0.4, k = 7.25, and two different waiting times. 
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FIG. 2. Parametric plot of the integratea res]fioffse func- 
tion M{tm + T, tw) and the correlation function C{tm + r, t^,) 
for k = 7.25. a) Tf = 0.4, Triangles: t» = 100. Crosses: 
tw = 1000. Circles: tw = 39810. The two straight lines have 
slopes -1.0 and -0.62. b) T/ = 0.3, tw = 1000. Circles: 
tw = 10000. The straight lines have slopes —1.0 and — 0..45. 
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